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Abstract 

The nanoscopic structure and the stationary propagation velocity of (l + l)-dimensional sohd-on- 
sohd interfaces in an Ising lattice-gas model, which are driven far from equilibrium by an applied 
force, such as a magnetic field or a difference in (electro) chemical potential, are studied by an 
analytic nonlinear-response approximation [P. A. Rikvold and M. Kolesik, J. Stat. Phys. 100, 377 
(2000)] together with kinetic Monte Carlo simulations. Here we consider the case that the system 
is coupled to a two-dimensional phonon bath. In the resulting dynamic [K. Saito, S. Takesue, 
and S. Miyashita, Phys. Rev. E 61, 2397 (2000); K. Park and M. A. Novotny, Comput. Phys. 
Commun. 147, 737 (2002)], transitions that conserve the system energy are forbidden, and the 
effects of the applied force and the interaction energies do not factorize (a so-called hard dynamic) . 
In full agreement with previous general theoretical results we find that the local interface width 
changes dramatically with the applied force. However, in contrast with other hard dynamics, this 
change is nonmonotonic in the driving force. Results are also obtained for the force-dependence 
and anisotropy of the interface velocity, which also show differences in good agreement with the 
theoretical expectations for the differences between soft and hard dynamics. However, significant 
differences between theory and simulation are found near two special values of the driving force, 
where certain transitions allowed by the solid-on-solid model become forbidden by the phonon- 
assisted dynamic. Our results show that different stochastic interface dynamics that all obey 
detailed balance and the same conservation laws nevertheless can lead to radically different interface 
responses to an applied force. Thus they represent a significant step toward providing a solid 
physical foundation for kinetic Monte Carlo simulations. 

PACS numbers: 68.35. Ct 75.60.Jk 68.43.Hn 05.10.Ln 
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I. INTRODUCTION 



Moving internal boundaries or interfaces separating different regions are present in many 
problems in nature, and the challenge of understanding the dynamics of such processes has 
become increasingly important. In recent years considerable efforts have been made toward 
understanding the large-scale structures of growing interfaces.^ In contrast, there has been 
little work related to the microscopic and nanoscopic scales. This is surprising since the 
nanoscopic interface structure plays a crucial role in important interface properties such as 
mobility and catalytic and chemical activity. Technologically, as the sizes of the smallest 
man-made structures decrease, interfacial properties become essential and even dominant. 
Nanoscale assemblies with highly ordered building blocks, such as quantum dots-'^ and 
quantum wires, must be fabricated on a surface or through an interface. 

The basic mechanisms of interface growth are complex and often unknown. A standard 
way to deal with this problem is constructing a stochastic model that reproduces essential 
features. However, extreme care has to be taken with this approach. Recent studies indicate 
that different stochastic dynamics, even when they have the same conserved quantities and 
satisfy detailed balance, lead to important differences in the nanostructure of field-driven 
interfaces.-i^'^i^i^ Surfaces driven by hard dynamics (in which the single-site transition rates 
cannot be factorized into one term that depends only on the interaction energies and a second 
term that depends only on the field energies, in contrast with soft dynamics for which 
this factorization is possible^), such as Glauber, Metropolis, and the two-step transition 
dynamics approximation (TDA),—^ have a strong dependence on the applied field. For all 
hard dynamics studied so far, the average step height increases dramatically with increasing 
field. In contrast, interfaces driven by soft dynamics, such as the soft Glauber- and the one- 
step-dynamics (OSD),-ii^ii^ are at most only weakly dependent on the field and relatively 
smooth. Furthermore, interfaces driven by hard dynamics, such as Glauber^ and TDA,- 
display significant asymmetry between the spin populations on their leading and trailing 
edges, while interfaces moving under soft dynamics either display no (soft Glauber-) or only 
weak (OSD^) anisotropy. 

In this paper we study by kinetic Monte Carlo (MC) simulation and a dynamic mean-field 
approximation the motion of a Burton-Cabrera-Frank solid-on-solid (SOS) interface^^ that 
evolves under a non- conservative dynamic resulting from coupling the system to a phonon 
heat bath. SOS interfaces belong to the Kardar-Parisi-Zhang (KPZ) dynamic universality 
class,"^'^^ in which the macroscopic, stationary distribution for moving interfaces is Gaus- 
sian, corresponding to a random walk with independent increments. The phonon-assisted 
dynamic is obtained by introducing a weak, linear coupling between a square-lattice Ising 
quantum ferromagnet and a phonon (i.e., bosonic) heat bath attached to the spin system. 
The transition rates have been calculated using the quantum-mechanical density matrix 
equationF^ii^ii2ii^i2i and most recently also by the lattice-frame method.-^^ Both methods 
give consistent results. The resulting dynamic is quite different from the Glauber dynamic, 
which can be similarly derived from coupling to fermionic baths.— In particular, for phonon 
baths of dimension greater than one, the phonon-assisted dynamic prohibits transitions that 
conserve the system energy, even if they are allowed by the SOS restriction. As a result, 
the model becomes non-ergodic near special values of the driving field, and the interfaces 
can get stuck in metastable states. The average step height and propagation velocity there- 
fore become nonmonotonic functions of the field. Phonon-assisted dynamics are relevant 
in a great variety of physical phenomena, ranging from the non-linear optical response of 
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semiconductors^^ to the dynamics of quantum dots.-i^ The derivation of a phonon-assisted 
stochastic dynamic is a significant step toward putting kinetic MC on a sohd physical foun- 
dation. 

In this paper we derive analytic, approximate expressions for the interface propagation 
velocity as a function of field, temperature, and interface orientation. Our approach is 
based on a mean-field approximation that assumes that individual steps on the interface 
are statistically independent ,->^i^ i.e., short-range correlations are neglected. This limitation 
becomes apparent when we compare the analytical and simulated results. 

The remainder of this paper is organized as follows. In Sec. [TTl we introduce the SOS 
interface model and give the transition rates for the phonon-assisted dynamic. Also in this 
section we summarize the mean-field approximation for the time evolution of the single- 
step probability density function (pdf), as well as its stationary form. We further give 
expressions for the spin-class populations and interface velocity in terms of the applied field, 
the temperature, and the angle of the interface relative to the lattice axes. In Sec. IIIII we 
compare simulations and analytical predictions for the detailed stationary interfacial nanos- 
tructure, including the asymmetry of the simulated nonequilibrium interfaces. A summary 
and conclusions are provided in Sec. [IVl 



II. MODEL AND DYNAMICS 

The SOS interfaces are described by the nearest-neighbor S = 1/2 Ising Hamiltonian 
with anisotropic, ferromagnetic interactions Jx and Jy in the x and y direction, respectively: 

= — ^ Sx,y {JxSx+l,y + JySx,y+l + H) , (l) 

x,y 

where Sx,y = ±1, Hx,y runs over all sites, and the applied field H is the driving force. 
The interface is introduced by fixing Sx^y = +1 and —1 for large negative and positive 
respectively. Without loss of generality we take H > 0, such that the interface on average 
moves in the positive y direction. This Ising model is equivalent to a lattice-gas model with 
local occupation variables Cx^y G {0, 1}. 25-26 Specifically, we identify s = +1 with c = 1 
(occupied or "solid") and s = —1 with c = (empty or "fluid"). 

The SOS model considers an interface in a lattice gas or S* = 1/2 Ising system on a square 
lattice of unit lattice constant as a single-valued integer function h{x) of the x-coordinate, 
with steps 6{x) = h{x + 1/2) — h{x — 1/2) at integer values of x. A typical SOS interface 
configuration is shown in Fig. [1] In this paper the two possible states of the site [x, y) are 
denoted by the two Ising spin values Sx^y = ±1. (In order that the step positions and the 
interface heights be integer as stated above, we place the spins at half-integer values of x 
and y, i.e., at the centers of the unit cells separated by dotted lines in Fig. [H) 

The interface is made to evolve under the phonon-assisted dynamic, a single-spin-fiip 
(nonconservative) set of transition rates that satisfy detailed balance for the allowed transi- 
tions. In most cases this ensures the approach to equilibrium, which in this case is a uniformly 
positive phase with the interface pushed off to positive infinity. (For exceptions, see below.) 
The dynamic is defined by the single-spin transition rates, W{sx,y ~^x,y ) = W{(3AE). 
Here /3 is the inverse of the temperature T (Boltzmann's constant is taken as unity), and AE 
is the energy change corresponding to a successful spin flip. The detailed-balance condition 
(valid for transitions between allowed states) is expressed as W{(3AE)/W{—PAE) = e~^^^ . 
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The transition rates for the phonon-assisted dynamic are defined a9^>i^ 

(2) 

where d G {1,2,3} is the dimension of the bosonic heat bath. Physically, this rate is the 
product of three factors: the phonon occupation number (e^''^^' — the phonon density of 
states, proportional to |A_E'|'^~^, and the magneto-elastic spin-phonon coupling, proportional 
to \/^E\M^ For = 2 and 3, W^b{T,/\E = 0) = 0, while for = 1 the transition rate 
is nonzero and smooth at AE = 0. These transition rates are plotted in Fig. [2](a). For 
comparison we also plot the transition rates for the Glauber dynamic in Fig. [2]^b). In the 
present work, with the exception of Fig. [7] we use d = 2, i.e, a two-dimensional heat bath. 
Thus, the transition rates vanish linearly with \AE\ near \AE\ = 0. It should, however, be 
emphasized that the derivations of Eq. ([2]) are based on a weak, linear coupling of the phonon 
bath to the spin system. It is therefore possible that nonlinear and/or multiphonon effects 
may set a lower bound on physical transition rates for AE near zero, and thus restore 
the ergodicity of the spin model. However, we note that recent experiments on phonon- 
mediated spin relaxation in a quantum dot shows a significant decrease in the relaxation 
rate for transitions involving AE near zero.'* 

Notice that the phonon-assisted transition rates cannot be factorized into one part that 
depends only on the interaction energy and another that depends only on the applied field; 
thus it belongs to the class of dynamics defined as hard.^'^- In order to preserve the SOS 
configuration at all times, fiips are allowed only at sites which have exactly one broken bond 
in the y direction. 

With the Ising Hamiltonian, there are only a finite number of different values of AE. The 
spins can therefore be divided into ten classes,— i^^i^Sii^ labeled by the spin value s and the 
number of broken bonds between the spin and its nearest neighbors in the x and y directions, 
j and k, respectively. The spin classes consistent with the SOS model are denoted jks with 
j e {0, 1, 2} and k e {0, 1}. They are shown in Fig. □ and listed in Table E At H = 0, 
AE = for transitions between 11— and 11-|- (diffusion of steps of unit height). Thus these 
transitions are forbidden for d = 2 and 3. At H = 23^, the transitions forbidden for d = 2 
and 3 are between 01— and 21-|- (nucleation or elimination of a knob of stable phase on 
a smooth, horizontal interface). For other values of (nonnegative) H, no transitions allowed 
by the SOS condition are forbidden. 

In the SOS model and our analytical approximation the heights of the individual steps are 
assumed to be statistically independent and identically distributed. This assumption is exact 
for H = 0.— The step-height probability density function (pdf) is given by the interaction 
energy corresponding to the \S{x)\ broken Ja,-bonds between spins in the columns centered 
at {x - 1/2) and {x + 1/2) as 

p[6{x)] = Z(0)-ixl'^(")l e^(<^)^(") . (3) 

The factor X determines the width of the pdf, and 7(0) is a Lagrange multiplier which 
maintains the mean step height at an x-independent value, {S{x)) = tan0, where is the 
overall angle between the interface and the x axis. Z{(f)) is a partition function that will 
be discussed below. In equilibrium, X is simply the Boltzmann factor, e~^^"^^, which is 
independent of H. In previous papers^i^ an expression for a field-dependent X(T,H) was 
obtained, based on a dynamic mean-field approximation for the equation of motion for the 
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single-step pdf together with a detailed-balance argument for the stationary state. This 
improved non-linear response approximation gives (see Ref.0 for details of the calculation), 



^' ^ I W[P{-2H-4J,)]+WmH-4J,)] / ' 

which is independent of 7(0). The dependence on the specific dynamic is evident here by 
the presence of the transition rates associated with the reversal of a single spin, W{l3AE). 
For H = 0, X is reduced to its equilibrium value, X(T,0) = e~^'^"^"=. For soft dynamics, 
where the field and the interaction terms factorize, the if-dependence in Eq. cancels 
out, while for hard dynamics X has a nontrivial dependence on H. Early results indicated 
that the SOS interfaces generated with the soft Glauber dynamic are indeed independent of 
H.— However interfaces generated with the OSD dynamic, which is also soft, show a weak 
dependence on the interface structure of the field.- 
The partition function for the interface is 

+00 1 

^ E A-I<le-.* ^ ,3,^-=-^, (5) 



where 7(0) is given by 



e7W 



1 + tan0 + 



X2)2tanV + 4X2 



1/2 



2X(l + tan( 



(6) 



(see details in Refs. liandQ). 

The mean spin-class populations, {n{jks)), are all obtained from the product of the 
independent pdfs for 6{x) and S{x+1). Symmetry of p[(5(x)] under the transformation 
(x, 0, (5) — s> {—x,—(j),—6) ensures that {n{jk—)) = {n{jk+)) for all j and k. Numerical 
results illustrating the breakdown of this up-down symmetry for large H are discussed in 
Sec. mil The general expressions for the class populations are given in the third column of 
Table [H details of the calculation can be found in Ref. 0. 

Whenever a spin flips from —1 to the corresponding column of the interface advances 
by one lattice constant in the y direction. Conversely, the column recedes by one lattice 
constant when a spin flips from -|-1 to —1. The corresponding energy changes are given in 
the second column in Table [B Since the spin-class populations on both sides of the interface 
are equal in this approximation, the contribution to the mean velocity in the y direction 
from sites in the classes jk— and jk+ becomes 

{vyijk)) = W i(3AEijk-)) - W iPAEijk+)) . (7) 

The mean propagation velocity perpendicular to the interface becomes 

{v^{T,H,(^)) = cos(PY.{^Uks)){vyUk)) , (8) 

where the sum runs over the classes included in Table [H It was shown in Ref. that Eq. dHJ 
reduces to the results for the single-step^'^'^'^ and the polynuclear growth^'^'^ models 
at low temperatures for large and small 0, respectively. The spin-class populations listed 
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in Table U can be calculated explicitly by replacing X with its corresponding value from 
Eq. Q. 

In the next Section we show that the nonlinear-response approximation gives good agree- 
ment with MC simulations of driven SOS interfaces evolving under the phonon-assisted 
dynamic for a wide range of fields and temperatures. The main deviations between theory 
and simulations are seen for H/ J = and 2, where some transitions allowed by the SOS 
restrictions have AE = and thus are forbidden by the two-dimensional phonon-assisted 
dynamic, Eq. ([2]). 

III. COMPARISON WITH MONTE CARLO SIMULATIONS 

We calculated the step-height distributions, propagation velocities, and spin-class popu- 
lations, analytically and by kinetic MC simulations, for the phonon-assisted dynamic in the 
isotropic case, Jx = Jy = J- The details of our particular implementation of the ra-fold way 
rejection-free MC algorithm^i^ are described in Refs. |^ andQ. The extension to continuous 
time, which is necessary to handle transition rates greater than unity, was introduced in 
Ref. i 

The numerical results presented here are based on MC simulations mostly at the two 
temperatures, T = 0.2Tc and 0.6Tc (Tc = —2 J/ ln(v^ — 1) ^ 2.269 J is the critical temper- 
ature for the isotropic, square- lattice Ising model^), with = 10 000 and fixed between 
and 45°. In order to ensure stationarity we ran the simulation for 10^ n-fold way updates 
per updatable spin (UPS) for thermalization before taking any measurements. Unless oth- 
erwise noted, the initial condition before thermalization was a microscopically flat interface. 
The initial condition only makes a difference near H = 0. Stationary class populations and 
interface velocities were averaged over 10^ UPS. For the stronger fields at T = 0.2Tc we used 
ten times as many UPS. Adequate statistics for the step-height pdfs were ensured by the 
large L^. 

A. Stationary single-step probability densities 

Stationary single-step pdfs were obtained by MC simulation at T = 0.2Tc and 0.6Tc for 
= and several values of H. The simulation data and the theoretical results for p[5] 
are shown in Fig. [3l The theoretical results are calculated with Eq. ([3]), with X{T, H) from 
Eq. (111). The agreement between theoretical and simulated results is quite good, particularly 
at the lower temperature. 

Another way to compare the analytical and simulation results is by calculating {\6\) by 
averaging over the simulated step-height pdf, and comparing these values with the theoretical 
ones obtained from Eq. ([3]), {\6\) = 2X/ (1 — X^), with X from Eq. (jl]). The results are 
shown in Fig. IHfor = at T = 0.2Tc and O.GTc, calculated theoretically (solid lines) and 
by MC simulation (symbols). The agreement between both results is reasonable. However, 
the theoretical data present a smoother dependence on the field than obtained from the 
simulation. 

In Fig.m the behaviors of {\6\) near H/ J = and 2 are of particular interest. At H/ J = 
the system should be in equilibrium, and the theoretical results are exact.— The discrepancy 
between theory and simulation at this field, especially at T = O.GTc, therefore means that 
the system simulated with the phonon-assisted dynamic, starting from a microscopically flat 
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initial state, does not equilibrate completely. This is not due to a too short thermalization 
time. Rather, the reason is the aforementioned suppression by Eq. ([2]) near H = of 
transitions between states 11+ and 11—, which correspond to diffusion of steps along the 
interface and represent an important mechanism for equilibration. However, when 
near H/ J = at T = 0.6Tc, is obtained by the phonon-assisted dynamic starting from 
the thermalized interface generated with the standard Glauber dynamic, there is excellent 
agreement with the theoretical result, as can also be seen in Fig. |H As H/J is increased 
from zero, the effect of the initial condition rapidly vanishes. 

In Fig.[5]^a) we show together snapshots of stationary interfaces at H/ J = for T = O.GTc, 
generated in three different ways: using the standard Glauber dynamic; with the phonon- 
assisted dynamic starting from the equilibrated interface obtained by the Glauber dynamic; 
and with the phonon-assisted dynamic starting from a microscopically flat interface. Due 
to the fact that energy-conserving moves (horizontal or vertical step-diffusion) are prohib- 
ited by the phonon-assisted dynamic, the interface started from the equilibrium interface 
is highly correlated with the latter, and both have {\6\) ~ 0.49, the equilibrium value. For 
the same reason, the phonon-assisted interface started from a microscopically fiat interface 
configuration does not fully equihbrate, but settles into a metastable configuration with 
^ 0.41, as seen in Fig. H 

At H/ J = 2, MC simulations give a value of {\S\) strictly zero for T = 0.2Tc and very close 
to zero for T = 0.6Tc, while the theoretical value is very small but nonzero at T = 0.2Tc, 
and clearly larger at T = O.GT^. (See inset in Fig. HI) For strong values of the field, the 
step height is only weakly dependent on the temperature. These results are quite different 
from those obtained with the standard (hard) Glauber dynamic (see Fig. 5(a) of Ref. 0). 
However, the strong H dependence of the step heights is characteristic of hard dynamics. 

In Fig. [5]^b) we show together snapshots of stationary interfaces at H/ J = 2 and T = 
O.GTc, one thermalized with the standard Glauber dynamic, and the other generated by the 
phonon-assisted dynamic, using the Glauber interface as its starting state. The interface 
obtained with the phonon-assisted dynamic is almost entirely microscopically fiat, with a 
very small density of "backward" (21—) notches that are created at a rate oc exp(— 8J/T) 
and annihilated almost immediately. In fact, for H/J = 2 the interface gets stuck as it 
can never progress beyond the absolute maximum of the starting configuration (and is thus 

non-ergodic) due to the vanishing rate of the transition 01 > 21+. One possible way to 

overcome this situation is to give the interface alternative paths to reach equilibrium. This 
could possibly be done by relaxing the SOS constraint to allow overhangs and bubbles.^-- In 
contrast, the Glauber interface at the same field and temperature propagates at a nonzero 
velocity and is microscopically quite rough, with {\6\) ~ 1.79. 

B. Stationary interface velocities 

In Fig. [6] we show the mean propagation velocity perpendicular to the interfaces vs H/ J 
for = 0, obtained with the analytical approximation, Eq. ([8]), and by simulations. In 
general there is good agreement between the MC results and the nonlinear-response theory. 
However, there is a significant discrepancy at H/ J = 2. At this field, the simulated velocities 
are zero, independent of the temperature, while the theoretical value is small but nonzero 
and increases with temperature. This is perfectly consistent with the microscopically fiat 
interface structure at H/J = 2, discussed above. In fact, the knob-nucleating transition 
01 > 21+, which is forbidden at this field, is precisely the transition needed to nucleate 
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the advance of a microscopically flat interface. When H/ J > 2, the velocity increases 
rather rapidly with H. This behavior is very different from the one obtained for other hard 
dynamics such as the standard Glauber and the TDA, where the velocity is bounded by 
unity. (See Fig. 6 of Ref. and Fig. 6(a) of Ref. i.) 

In contrast to the results discussed above, Fig. [7] shows the velocity for the case in which 
the dimension of the phonon bath is unity, i.e., c? = 1, which has a transition rate that 
decreases smoothly and monotonically with /S.E/T (in contrast with the d = 2 and 3 cases 
in which the transition rates vanish for H = and H/J = 2). The agreement between 
theory and simulation is excellent over the whole range of H/ J and tilt angles, except for 
very large angles at higher fields. (This is also the case for step-height distributions and 
other characteristic quantities.) 

The dependence of the normal velocity on the tilt angle (p is shown in Fig. [8]^a) and 
Fig. El^b) for several values of H/ J aXT = 0.2Tc and T = 0.6Tc respectively. The agreement 
between the theoretical results and the simulations is very good except at higher fields, 
where the agreement is only good at intermediate values of 0. The results are qualitatively 
similar to those obtained with other hard dynamics (see Refs.0and §). At T = 0.2Tc, in 
weak fields the velocity increases with 0, in agreement with the polynuclear growth model 
at small angles and the single-step model for larger angles. For strong fields the behavior 
changes gradually to the reverse anisotropy of Eden-type models.—*^ This is essentially 
the same behavior observed for the TDA dynamic- At T = 0.6Tc, the velocity is nearly 
isotropic for weaker fields, while becoming Eden like for stronger fields. The exception is 
the case H/J = 2, which at small angles presents a polynuclear-growth type, as well as 
significant differences between the theoretical and simulated results. The behavior of the 
normal velocity at T = 0.6Tc (excluding the case H/ J = 2), is very similar to that observed 
for both the TDA and the OSD dynamics.- 

The temperature dependence of the normal interface velocity is shown in Fig.[9]for several 
values of H/ J. The agreement between the simulations and the analytical results is reason- 
able, except for H/J = 2, where the simulated velocity remains zero for all temperatures, 
while the predicted velocity increases monotonically with the temperature. This discrep- 
ancy is also due to the metastable, static and microscopically fiat, interface that forms at 
H/J = 2. This figure also shows that as T 0, the system develops a step discontinuity: 
the velocity is zero for H/ J <2 and increases with H for stronger fields. This discontinuity 
at T = is also observed with the TDA^ and with the standard Glauber dynamic (see Fig. 8 
of Ref. 0). 

C. Spin-class populations and skewness 

To test the analytical assumption that different steps are statistically independent, we 
compare the analytical results for the mean class populations^ with the simulated ones. The 
six mean class populations — (n(Ols)), (n(lls)), and (n(21s)) with s = ±1 — for = at 
T = 0.2Tc and T = 0.6Tc are shown vs H in Fig. [TOl At both temperatures, the analytical 
approximations follow the average of the populations for s = +1 and s = — 1 qualitatively 
well. However, for small fields, H/ J < 2, the simulations show a stronger dependence on 
H than the mean-field results. This difference is more evident at the higher temperature, 
where the simulations show that the population in front of the interface (s = —1) is quite 
different from the one behind it (s = +1). Well away from the special fields {H/J = 
and 2), the interfaces are a little rougher than the theory predicts (lower 01s and higher 
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lis populations). Near the special fields the interfaces appear to get caught in smoother 
(metastable) configurations. For H = 0, the 01s populations are slightly higher that the 
theoretically predicted value, and the lis populations are slightly lower. This is more 
evident at the higher temperature. Note that when H/ J = the interface should be in 
equilibrium, and the theoretical value is exact. However, this exact value is only reached 
if the interface is equilibrating properly. This can be seen in Fig. [TUT b) where we include 
the population averages, both calculated starting from a microscopically flat interface and 
from a thermalized interface generated with the standard Glauber dynamic. For H/ J = 2 
the interface is much smoother than predicted. (For T = 0.6Tc, the measured values for 
(n(01±)), (r;,(ll— )), and (r;,(ll+)) are approximately 0.9999, 10~^, and 0, compared with 
the respective predicted values of 0.9, 0.1 and 0.01.) 

The short-range correlations between neighboring steps are responsible for the skewness 
between the spin populations on the leading and trailing edges of the interface that appears 
in the simulation results. This phenomenon is commonly observed in driven interfaces. It 
occurs even when the long-range correlations vanish as they do for interfaces in the KPZ 
dynamic universality class, to which the present model belongs for all finite, nonzero values 
of H/ J 7^ 2. Skewness has also been observed in several other SOS-type models, such 
as the body-centered SOS model studied by Neergaard and den Nijs,— the model for step 
propagation on crystal surfaces with a kink-Ehrlich-Schwoebel barrier studied by Pierre- 
Louis et al.,-^ and a model for the local time horizon in parallel kinetic MC simulations 
studied by Korniss et al.— No skewness was observed for the SOS model with the soft 
Glauber dynamic- However, some skewness was present in the OSD model,- indicating that 
a complete lack of skewness is not a necessary characteristic of soft dynamics. Also, a small 
degree of skewness was observed for the Ising model (whose interfaces include bubbles and 
overhangs) with the soft Glauber dynamic (about two orders of magnitude smaller than 
the skewness observed for the hard Glauber dynamic).^ The TDA dynamic also presents 
considerable skewness.^ The correlations associated with the skewness generally lead to 
a broadening of protrusions on the leading edge ("hilltops"), while those on the trailing 
edge ( "valley bottoms" ) are sharpened,— or the other way around.— In terms of spin-class 
populations, the former corresponds to (n(21— )) > (?t,(2H-)) and (n(lH-)) > (n(ll— )). 
The relative skewness can therefore be quantified by the two functions, '^^ 

_ (n(21-)) - (n(21+)) 
^ M21-)) + (n(21+)) ' 

and^ 

(n(ll+)) - (n(ll-)) 

(n(ll+)) + (n(ll-)) • ^ ^ 

These two skewness parameters are shown together in Fig. [TTl The temperature dependence 
of the skewness is stronger at the lower temperature and smaller fields, and it is especially 
pronounced for p, due to the low concentration of sites in the class 21+ at low temperatures 
and weak fields, e is very small and almost independent of T and H, except at H/J = 2.0, 
where both values are near unity, consistent with the picture of an interface with a very low 
density of 21— notches. 
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IV. DISCUSSION AND CONCLUSIONS 



In this paper we have studied the nanostructure of an unrestricted SOS interface inter- 
acting with a two-dimensional phonon heat bath and driven far from equihbrium by an 
apphed field. This work is a continuation of previous studies aimed to explore the crucial 
role of the stochastic dynamics selected to simulate physical systems.-'^!^'^'^'^'^ Important 
properties such as nanostructure and mobility of driven interfaces have been shown to be 
strongly dependent on the stochastic dynamics.-'^'^'^'^ 

The transition rates that give the evolution of the system are calculated by coupling 
the spin Hamiltonian to a two-dimensional phonon heat bath. The dynamic generated in 
this way belongs to the class known as hard.— In condensed-matter systems, phonon-driven 
dynamics are dominant, and technologically they are becoming increasingly important, e.g., 
in connection with the behavior of quantum dots.- 

We studied the nanostructure and velocity of the SOS interface by kinetic MC simulations 
and by a non-linear mean-field theory developed in previous papers.™*^ We calculated the 
interface velocity as a function of the driving field, temperature, and angle of the interface 
relative to the lattice axes. We also studied the local shape of the interface in terms of the 
spin-class populations, the average height of a step, and the probability density for individual 
steps in the interface. 

In general we found good agreement between the theoretical calculations and the MC 
simulations. In particular we found the strong dependence of the interface structure on the 
field characteristic of systems that evolve under hard dynamics, such as the Glauber or TDA 
dynamics.-i^i^ 

Our theoretical results are based on the mean-field assumption that individual steps of 
the interface are statistically independent; short-range correlations are neglected. However, 
our MC results show asymmetry between the spin populations on the leading and trailing 
edges of the interface, which is an indication of the existence of such short-range correlations. 
With increasing field, the interfaces undergo a gradual breakdown of up-down symmetry, 
which has also been observed in other examples of driven interfaces.-i2^>^i^ Aside from 
such, relatively minor, discrepancies between the theoretical mean-field predictions and the 
simulation results, which show that there is room for improvement of the mean-field model, 
the theory predicts very accurately the quahtative behavior of the interfaces and yields a 
reasonable over-all approximation to their quantitative behavior. The important exceptions 
are the special field values, H/J = and 2, where certain transitions allowed by the SOS 
constraint are forbidden by the phonon-assisted dynamic for phonon baths of dimension 
greater than one. At H/ J = this leads to a failure of the simulated interface to reach 
thermal equilibrium, while at H/ J = 2 it leads to a metastable, abnormally flat interface 
that is unable to propagate. 

It should, however, be noted that the phonon-assisted dynamic defined by Eq. ([2]) is 
based on a weak, linear coupling of the bosonic bath to the spins. It is therefore possible 
that higher-order and/or multiphonon corrections to the transition rates could restore the 
vanishing rates for energy-conserving transitions. Nevertheless, higher-order effects may not 
completely mask the slowing-down of the interface in the vicinity of the special field values 
observed here. We therefore expect that much of the characteristic field dependence will 
carry over to more sophisticated rate models. This expectation is supported by the results 
of recent experiments on phonon- mediated spin dynamics in a quantum dot,- in which 
significant slowing-down was observed for nearly energy- conserving transitions. Another 
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interesting question is to what extent a relaxation of the SOS constraint to consider a full 
Ising interface including overhangs and bubbles (see Refs. i andlsl) might open alternative 
channels for full equilibration. Naturally, one can also think of more general dynamics, 
including multiple-spin-flip elementary transitions. These questions are left for future study. 

As in previous studies our results indicate strong differences between inter- 
faces moving under different stochastic dynamics, emphasizing the need for extreme care in 
selecting the appropriate dynamic for the physical system of interest. This general under- 
standing and the specific results for the phonon-assisted dynamic presented in this paper 
represent significant steps in the direction of putting kinetic MC simulations on a solid 
physical foundation. 
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TABLE I: The spin classes in the anisotropic square-lattice SOS model. The first column contains 
the class labels, jks. There are two other classes, 10s and 20s, that also have nonzero populations 
in the SOS model, but are not included because flipping a spin in any of them would produce an 
overhang or a bubble and is therefore forbidden. The second column contains the change in the 
total system energy resulting from reversal of a spin from s to — s, AE(jks). The third column 
contains the mean spin-class populations for general tilt angle (p, with cosh7(i;^) from Eq. 



Class, jks 


AE{jks) 


(n(jfcs)) 


01s 


2sH + 4 Jj; 


1-2X cosh7(fli)+X^ 


(1-X^)^ 


lis 


2sH 


2X[(l+X^)cosh7((^)-2X] 


(1-X2)2 


21s 


2sH - 4 J^. 


X'-' [1-2X cosh 7(fli)+X2] 


(1-X^)^ 



12 













00- 






21+ 








01- 




20+ 




11- 


21+ 


11+ 1 21- 


20+ 


20- 


11+ 


10+ 




01+ 


10+ 


21- 


10+ 




00+ 






01 + 







■^-3 -2-10 12 3 



X 

FIG. 1: A short segment of an SOS interface y = h{x) between a positively magnetized phase 
(or "sohd" phase in the lattice-gas picture) below and a negative (or "fluid") phase above. The 
step heights are 8{x) = h{x + 1/2) — h{x — 1/2). Interface sites representative of the different SOS 
spin classes (see Table U ) are marked with the notation jks explained in the text. Sites in the 
uniform bulk phases are 00— and 00+. This interface was generated with a symmetric step-height 
distribution, corresponding to = 0. 




FIG. 2: (a) (Color online) The transition rates for the d-dimensional phonon-assisted dynamic, 
WpB, shown scaled by T'^ vs the energy difference A£^ scaled by T. (b) The transition rates for 
the standard Glauber dynamic, Wq,, shown vs /S.E/T. 
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FIG. 3: (Color online) MC (data points) and analytical (solid lines) results for the stationary 
single-step pdf, shown on a logarithmic scale vs <5, for the values of H/ J given in the legend, (a) 
T = Q.2Tc. (b) T = 0.6Tc. The symbols (and colors) have the same interpretations in (a) and (b). 
Note the nonmonotonic field dependence near H/ J = 2. 




FIG. 4: (Color online) Average stationary step height {\6\) vs H/ J for 0=0 at T=0.2Tc and 0.6Tc. 
The curves represent the theoretical results. Curve with circles (black): T = 0.2Tc. Curve with 
squares (gray, red online): T = 0.6Tc. These data refer to interfaces started from a microscopically 
flat initial state. We also include values near H/ J = 0, calculated by the phonon-assisted dynamic 
using as starting state the thermalized interface obtained with the standard Glauber dynamic. 
Black asterisks: T = 0.2Tc. Grey asterisks (red online): T = O.GTc. The differences are only 
evident near H/J = 0. In this and all the following figures, the statistical uncertainty is much 
smaller than the symbol size. The inset shows a magnified view of the region around H/J = 2. Note 
the disagreement between the theoretical and the simulation values at T = 0.6Tc when H/J = 2. 
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FIG. 5: (Color online) Short segments of thermalized interfaces, (a) H/J = 0, T = O.GTc. 
The three graphs show an equilibrium interface created with the standard Glauber dynamic over 
10^ UPS (medium gray, red online); an interface created by the phonon-assisted dynamic over 
10® UPS, using the equilibrium interface as starting state (dark gray, blue online); and an interface 
created by the phonon-assisted dynamic over 10^'^ UPS, using a microscopically flat interface as 
starting state (black), (b) H/J = 2, T = 0.6Tc. The jagged interface (medium gray, red online) is 
in the statistically stationary state, propagating in the direction of the arrow under the standard 
Glauber dynamic. At a given time, the dynamic is switched to the phonon-assisted transition rates, 
using the Glauber interface as initial state. The lagging parts rapidly catch up with the absolute 
maximum of the Glauber interface, where the interface gets permanently stuck in an almost perfect, 
microscopically flat configuration with a very small density of "backward" 12— notches (circled). 
At this field, the transition forbidden by the phonon-assisted dynamic is the nucleation of "forward" 
12+ notches, which are needed to nucleate propagation of a microscopically flat interface. 
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FIG. 6: (Color online) The average stationary normal interface velocity {v±) vs H/J for (f) = 0. 
The MC results are shown as data points, circles for T = 0.2Tc and squares for T = 0.6Tc, and the 
theoretical results as solid curves. The inset shows a magnified view of the region around H/J = 2. 
Note again the disagreement between the theoretical and the simulation values at T = 0.6Tc when 
H/J = 2. 
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FIG. 7: (Color online) Average stationary normal interface velocity at T = 0.2Tc, obtained by 
coupling the system to a one- dimensional phonon bath. MC data are represented by the symbols 
and analytical results by the solid curves, (a) Velocity vs H/J for (j) = 0. (b) Velocity vs tan (p for 
several values of H/J. 




FIG. 8: (Color online) The average stationary normal interface velocity {v±) vs tan0, for several 
values of H/J. The symbols represent MC data, and the solid curves analytical results, (a) 
T = 0.2Tc, (b) T = 0.6Tc. The symbols have the same interpretations in (a) and (b), given by the 
legend in (a). Online, the colors of the curves and symbols match. 
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FIG. 9: (Color online) The average stationary normal interface velocity {v±) vs T for = and 
H/J between 0.5 and 3. MC data are represented by data points and analytical results by solid 
curves. Online, the colors of the curves and symbols match. The agreement between simulation 
and theory is quite good, except for H/J = 2. 




FIG. 10: (Color online) Mean stationary class populations {n{jks)) vs H/J for = 0. The 
simulation results are indicated by symbols, and the analytic approximations by solid curves, (a) 
T = 0.2Tc. (b) T = O.GTc. At T = 0.6Tc we also include some values (asterisks for {n{jk+)) and 
pluses for {n{jk—)) ) calculated with the interface created by the phonon-assisted dynamic, using 
the thermalized interface obtained by the Glauber dynamic as starting state. Note that in this 
case there is excellent agreement between theory and simulations at H/J = 0. The other symbols 
have the same interpretations in (a) and (b), given by the legend in (a). The insets in both (a) 
and (b) show (n(lls)) and (n(21s)) near H/J = 2. 
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FIG. 11: (Color online) The two relative skewness parameters p (circles, black) and e (squares, red 
online), defined in Eqs. ([9|) and pop . respectively. The parameters are shown vs H/J for = 0, at 
T = 0.2Tc (empty symbols) and at T = O.GTc (filled symbols). The inset shows a magnified view 
of the region around H/J = 2. 
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